### This script will produce a BRT model of underlog temperature
### Model parameters were determined by an optimization procedure
		
# Load Elith Source Code

setwd('/home1/99/jc152199/brt/')
source('brt.functions.R.cjsedit.r')

### Define parameters

tc = ''
tf = 1
nt = 7000
lr = ''

# Load gbm library

library('gbm')

#### Reset work directory

output.dir = '/home1/99/jc152199/brt/underlog/emptempvarsandlogvars/'

setwd(output.dir)

# Read in data to model, which consists of the training and testing set of the optimal model

train.data = read.csv('/home1/99/jc152199/brt/underlog/emptempvarsandlogvars/p_02/UL_BRT_Model_Data_Plus_Preds.csv',header=T)
test.data = read.csv('/home1/99/jc152199/brt/underlog/emptempvarsandlogvars/p_02/UL_BRT_Test_Data_Plus_Preds.csv',header=T)

### Rowbind train.data and test.data to recreate model.data

model.data = rbind(train.data,test.data)

### Remove some objects

rm(train.data)
rm(test.data)

#### Remove site-days with residual values greater than....
#### First calculate residuals

model.data$resid = model.data$ulmax - model.data$ulmaxpreds

#### Find residuals greater than 3.75 (only want to remove 1% of worst residual values)

model.data = model.data[which(abs(model.data$resid)<3.75),]

# Randomly shuffle model.data before subsetting to training/testing sets

model.data = model.data[c(sample(c(1:nrow(model.data)), nrow(model.data),replace=FALSE)),]

#### Remove sites with both a log density and circumference

model.data = model.data[which(is.na(model.data$logden)==F & is.na(model.data$logcirc)==F),]

# Randomly sample half of data set to produce a training and a testing set

sample.vector = c(sample(c(1:nrow(model.data)), round(nrow(model.data)*as.numeric(tf),0),replace=FALSE))

### Create training and testing sets

train.data = model.data[c(sample.vector),]
test.data = model.data[-c(sample.vector),]

# Run BRT model using gbm.step, parameters are defined as arguments, indepedent variables are limited to BRT max, min, and range
		
brt.gbm.step = gbm.step(data=train.data,gbm.x = c(14,16:20),gbm.y = 3,max.trees=as.numeric(nt),family = "gaussian",tree.complexity = as.numeric(tc), learning.rate = as.numeric(lr) ,bag.fraction = 0.5, n.folds=10)

###  Save model results

save(brt.gbm.step, file=paste(output.dir,'ULModel.Rdata',sep=''))

# Produce partial dependence plots

gbm.plot(brt.gbm.step,variable.no = 0, nt = brt.gbm.step$n.trees, plot.layout = c(2,3))

# Append preds to model data

train.data$ulmaxpreds = brt.gbm.step$fitted

### Write training data (with preds) and testing data (with preds) to the output directory

write.csv(train.data, file=paste(output.dir,'UL_BRT_Model_Data_Plus_Preds.csv',sep=''),row.names=F)

# Write out model parameters and total deviance as a .csv file

params = data.frame(learning.rate = lr, tree.complexity = tc, train.fraction = tf, n.trees = nt, op.tree.num = brt.gbm.step$n.trees, train.deviance =  brt.gbm.step$cv.statistics$deviance.mean)

write.csv(params, file=paste('ul.maxtemp.model.summary.csv',sep=''),row.names=F)

### Make some plots from training data
#### Start with obs versus preds

### Start .png driver

png(paste(output.dir,'Training_Data_Obs_Versus_Preds.png',sep=''),units='cm',height=12, width=12, res=1000)

#### Plotting limits

lims = range(c(train.data[,3],train.data[,21]),na.rm=T)

#### Plot up testing obs versus preds

plot(train.data[,3],train.data[,21], xlab = 'Obs Training Data', ylab = 'Training Fitted Values', main = paste('Obs vs Preds - Deviance - ',round(brt.gbm.step$cv.statistics$deviance.mean,4),sep=''), xlim=lims, ylim=lims, col='olivedrab1')

#### 1:1 line

abline(a=0,b=1, lty=2, lwd=1, col='red')

### LM

lm1 = lm(train.data[,21]~train.data[,3]) # Perform a linear model using the Empirical Data as y-values and the preds (AWAP or microCLIM) as x-values

#### Plot lm relationship

abline(lm1,col='lightseagreen',lty=2, lwd=1)

#### Legend

legend('topleft',legend = c(paste('Adj. r^2 ',substr(summary.lm(lm1)[9],1,4),sep=''),paste('Slope ',round(lm1$coefficients[2],2),sep=''),paste('Intercept ',round(lm1$coefficients[1],2),sep='')),text.col=c('lightseagreen','lightseagreen','lightseagreen'), bty='n')

#### Shut device

dev.off()

#### Now plot obs versus resid

### Start .png driver

png(paste(output.dir,'Training_Data_Obs_Versus_Resids.png',sep=''),units='cm',height=12, width=12, res=1000)

#### Plotting limits

xlims = range(train.data[,3],na.rm=T)
ylims = range(train.data[,3] - train.data[,21],na.rm=T)

#### Plot up testing obs versus preds

plot(train.data[,3],train.data[,3]-train.data[,21], xlab = 'Obs Training Data', ylab = 'Training Resid Values', main = paste('Obs vs Preds - Deviance - ',round(brt.gbm.step$cv.statistics$deviance.mean,4),sep=''), xlim=xlims, ylim=ylims, col='olivedrab1')

### LM

lm1 = lm(train.data[,3]-train.data[,21]~train.data[,3]) # Perform a linear model using the Empirical Data as y-values and the preds (AWAP or microCLIM) as x-values

#### Plot lm relationship

abline(lm1,col='lightseagreen',lty=2, lwd=1)

#### Legend

legend('topleft',legend = c(paste('Adj. r^2 ',substr(summary.lm(lm1)[9],1,4),sep=''),paste('Slope ',round(lm1$coefficients[2],2),sep=''),paste('Intercept ',round(lm1$coefficients[1],2),sep='')),text.col=c('lightseagreen','lightseagreen','lightseagreen'), bty='n')

#### Shut device

dev.off()

### Done